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A small-scale turbulence model 

By T. S, Lundgren 1 


A model for the small-scale structure of turbulence is reformulated in such a way 
that it may be conveniently computed. The model is an ensemble of randomly ori- 
ented structured two dimensional vortices stretched by an axially symmetric strain 
flow. The energy spectrum of the resulting flow may be expressed as a time integral 
involving only the enstrophy spectrum of the time evolving two-dimensional cross 
section flow, which may be obtained numerically. Examples are given in which a 
Jfc- 5 /3 spectrum is obtained by this method without using large wavenumber asymp- 
totic analysis. The Jfc“ 5 / 3 inertial range spectrum is shown to be related to the 
existence of a self-similar enstrophy preserving range in the two-dimensional en- 
strophy spectrum. The results are insensitive to time dependence of the strain-rate, 
including even intermittent on-or-off strains. 


1. Introduction 

One of the issues in turbulence theory is to understand what kinds of elemen- 
tary flow structures are responsible for the part of the turbulent energy spectrum 
described by Kolmogorov’s celebrated law. A number of years ago, Lund- 

gren (1982) proposed a model of the small scale structure of turbulence which gives 
this spectrum. The model was put together as an ensemble of randomly oriented 
vortices with spiral structure, each vortex being subjected to an axially symmetric 
irrotational straining field. The strain was made constant and axially symmetric 
for analytical reasons and this should be thought of as a representation of much 
more complicated strains. The idea was to model the most important property of 
turbulence, namely the mixing property which causes fluid particles to rapidly sepa- 
rate, stretching vortex blobs into elongated tubes. This model generalizes an earlier 
model by Townsend (1951) which assumed randomly oriented Burgers vortices and 
gave a k~ l spectrum. 

Differential rotation in the vortices (the inner part has higher angular velocity 
therefore winds faster) causes the vortex layers in the spiral to tighten, and the axial 
straining decreases the cross section of the structure. These mechanisms cause a 
cascade to smaller scales which differs from the traditional concept of a cascade 
through eddies of different sizes induced by instabilities. 

Many turbulent flows have distinct two dimensional vortices, and often flow vi- 
sualization by laser-induced-fluorescence or smoke shows vortex cross sections with 
some spiral structure. Such structures may be seen in sections 4 and 6 of Van 
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Dyke’s (1982) book. Schwarz (1990) identified intermittently occurring layered vor- 
tex sheets (apparently spirals) in oscillating grid turbulence by flow visualization 
with small suspended crystalline platelets. On the other hand, the numerically 
simulated periodic box turbulence of Vincent and Meneguzzi (1991) contains very 
pronounced two-dimensioned vortices in which very little internal structure can be 
seen. (However, there is some spiral structure in their figure 15.) Perhaps the 
Reynolds number is too small in this kind of flow. 

A number of authors have made further studies based on the original paper by 
Lundgren. In particular, Lundgren (1985) applied the model to the calculation 
of the scalar spectrum of the product of a fast chemical reaction. Gilbert (1988) 
used similar ideas for the study of two-dimensional turbulence. Buntine and Pullin 
(1988) and Pullin and Buntine (1989) used the model as a basis for computations 
of spectra produced by the merger of vortices. Pullin and SafFman (1992) used it 
to calculate vorticity and velocity derivative moments for homogeneous isotropic 
turbulence. Recently, Gilbert (1992) has produced a qualitative cascade argument 
based on the model to explain the Kolmogorov spectrum. 

In section 2, the model is reformulated in such a way that the two-dimensional 
part of the model, the spiral flow or some more general two-dimensional vortical 
flow, may be conveniently carried out numerically. 

In section 3 numerical computations of two-dimensional flows are used to generate 
three-dimensional energy spectra. 

2. Reformulation of the spiral vortex model 

It will be useful to separate the strictly two-dimensional part of the model, the 
flow in the cross section of the vortex, from the straining part of the model which 
is performed by a transformation. The 1982 paper will be referred to as L and 
equations from that paper will be referred to by L:( ). It was shown in L that 
if a two-dimensional flow with vorticity given by u? 2 (x,y,tf) is placed in an axially 
symmetric strain flow with velocity components (— .5ax, —.5 ay,az) where the strain- 
rate a may be a function of time, then the vorticity in the resulting three-dimensional 
flow (which has the same initial vorticity as the two-dimensional flow) is given by 


«(*,V,0 = S(t ) o>2 (S(0 I/2 *,S(0 I/2 »,T(0) 

where 

/mi \ 

(2.1) 

S(t) — exp 

is the amount the flow is stretched and 

U 

(2.2) 

T(t) = 

[ S(t')dt' 

'o 

(2.3) 


is a strained time. The energy spectrum of an ensemble of strained vortices of all 
ages, accounting for the greater length of the older vortices, may be deduced from 
L:(58), L:(75) and L:(76) and expressed, in a new form, by 


E(k) = ^ j ' S(T) 1/2 F 2 (s(T)- 1/2 Jfc,r) dT 


(2.4) 
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where C = 2tt 2 l 0 N c /L 3 is a constant. The function F 2 is the enstrophy spectrum 
(the vorticity power spectrum) of the strictly two-dimensional flow. This is defined 


by 


where 


fin 

F 2 (k,t) = k I \ui 2 (k cos 9 k,k sin 8 k-,t)\ 2 ddk 

Jo 

w 2 = J j exp [~i(k z x + k 9 y)]u 2 (x,y,t)dxdy 


(2.5) 

( 2 . 6 ) 


and represents the enstrophy in a circular shell in wavenumber space, divided by 
the width of the shell. The integration variable in Eq.(2.4) is the strained time 
and the stretching function 5 must be expressed as a function of this strained 
time. When the strain-rate is constant, the usual case considered, S = exp(af) and 
T = (exp(af) — l)/a and hence 

5=1 + aT (2.7) 


is the appropriate function. In appendix A, it is shown by example that, even when 
the strain-rate is variable and quite different from constant, the stretching function 
is roughly a linear function of the strained time and the energy spectrum computed 
from Eq.(2.4) is insensitive to these modest deviations from linearity. 

In Eq.(2.4), the function F 2 is all that is needed from the two-dimensional flow, 
and it may be specified numerically or analytically. The time integral represents the 
effect of stretching. The finite time cut-off on the integral is to prevent the vortices 
from being stretched indefinitely. It was assumed in L that they ultimately coalesce 
into shorter vortices and renew the spiral structure. 

In L an analytical spiral vortex solution was developed. An approximate solution 
of the two-dimensional Navier-Stokes equation was given in the form of a Fourier 
series in the cylindrical angle variable 9, 


u(r,8,t) = ^2 u n (r,t)exp(in8) (2.8) 


w„(r,t) = /„(r)exp 

where the functions /„ are arbitrary and the angular velocity Q(r) is related to the 
average vorticity u>o by 

ruj 0 -- — r 2 fl(r). (2.10) 

dr 

The function fl(r) must be monotone decreasing - the property which gives differ- 
ential rotation. If the functions /„ are all the same (independent of n) the solution 
looks like a spiraling vortex sheet in the inviscid limit. The approximations in 
Eq.(2.8) require that t must be large, i.e., the error becomes small as t -+ oc. 
However, it will be seen numerically that it is quite good even for fairly small t. 


—infl(r)t 


1 2 /'dQY 

~ 3 n V d7 / 


vr 


(2.9) 
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By using this analytical solution, the enstrophy function F 2 may be expressed as 
a series of integrals which involve Bessel functions. These may be evaluated by the 
method of stationary phase, which requires both k and t to be large, an asymptotic 
result which may be written 


= t l G(k/t)exp(—2i/k 2 t/3) (2.11) 

where 

= < 212 > 

with 

k + nfl'(r n )t = 0. (2.13) 

The last equation results from the method of stationary phase. One is supposed 
to solve this for r n and substitute it into Eq.(2.12). Note that r n is a function 
of k/t , hence the form of the argument of the function G. In writing this result, 
the contribution from the n = 0 term has been omitted. This term has a different 
functional dependence and was shown in L to contribute little at high wavenumber. 
Equation (2.11) does not appear in L but may be deduced from L:(66) and L:(64) 
(with 5=1 and P = t). 

It is the functional form of Eq.(2.11) which is important. First note that the 
functions /„ are assumed to be of limited extent so that the spiral is restricted to 
a halo around a central vortex. Then the function F 2 looks like a localized hump 
when plotted versus wavenumber (this will be clear when some computations are 
seen). If the viscous factor can be neglected, the area under the hump stays constant 
because 

I t l G(kft)dk = constant. (214) 

Jo 

The similarity form of the function shows that the top of the hump moves to higher 
wavenumber with constant speed, while the width broadens and the peak decreases 
in such a way that the area stays constant. This is an enstrophy preserving temporal 
cascade. The physics is clear: as the spiral turns tighten due to differential rotation, 
the enstrophy shifts to higher wavenumber while being conserved. 

The function F 2 given by Eq.(2.11) is only part of the enstrophy spectrum. There 
is an additional large part at low wave number (the n = 0 term which was omitted). 
If viscosity is neglected, the total enstrophy is preserved, i.e., f u?|d A = J = 
constant. What has been shown here is that the spiral solution has scale separation 
of the enstrophy spectrum into two distinct peaks, and the enstrophy of each of 
these parts is independently preserved. 

The energy spectrum expression given by Eq.(2.4) could also have been cast in 
terms of the two-dimensional energy spectrum, E 2 say. For the spiral solution, £2 
also has scale separation into a large low wavenumber peak and a secondary smaller 
(much smaller because of the k~ 2 factor) high wave number peak with a similarity 
structure. While the total energy is conserved (in the absence of viscosity), the 
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energy in the separate peaks is not. The energy in the high wavenumber peak 
decreases like t“ 2 , giving up this energy to the low wavenumber peak. 

This enstrophy preserving similarity form is responsible for the A; -5 / 3 * part of the 
three-dimensional energy spectrum. When Eq.(2.11) is substituted into Eq.(2.4), 
the result 

E(k) = exp(— 2vk 2 /Zd) (2.15) 

is obtained if one takes the upper limit to be infinite and approximates Eq.(2.7) 
by 5 ~ aT, i.e., most of the contribution to the integral comes from large values 
of 5. Equation (2.15) results from a simple change of the integration variable to 
T/fc 2 / 3 and doesn’t depend at all on the specific form of G. This generalizes the 
model since there could be other kinds of flows with scale-separated self-similar 
conservative enstrophy spectra, although none are presently known. 

The following qualitative derivation of the functional form of the enstrophy spec- 
trum was motivated by a discussion with Javier Jimenez. Assume that the enstrophy 
spectrum has a self-similar form t _ft G(fc/^), then it is easy to see that conserva- 
tion of enstrophy implies a — ft. One may argue that in a fairly steady shear flow 
the length of an element of a vortex sheet increases linearly with <; its thickness, 
therefore, decreases like t~ l . The largest wavenumber in a system of such vortex 
sheets behaves like the reciprocal of the thickness. Therefore, k ~ t and a = f3 = 1, 
as desired. 

It can also be noted that if there is no axial straining, so that 5=1, Eq.(2.4) 
and Eq.(2.11) (without the viscous part) gives 

E(k) ~ i T 2 . (2.16) 

This was also noted by Gilbert (1988) and is consistent with Townsend’s (1951) 
observation that a fc“ 2 spectrum results from random arrays of vortex sheets. The 
conclusion to be drawn is that the axial straining gives the extra k 1 ^ 3 . 

The integration in Eq.(2.4) may be understood as follows. Imagine a A", T plane 
with S l ! 2 F 2 (K, T) plotted in a third dimension. This function looks like a ridge 
centered along a ray from the origin of the A', T plane, decreasing in height and 
spreading for increasing T. The argument K = 5“ 1 / 2 fc required in Eq.(2.4) is a 
curve in the A", T plane starting from K = k at T = 0 and moving to smaller I\ for 
larger T. The integrand in Eq,(2.4) is the height of the ridge seen as one traverses 
this curve. In order to get fc~ 5 / 3 for a band of fc-values, the integration must extend 
through the ridge into the small values on the far side, for all fc-values in this band. 

3, Computations 

The formula expressed by Eq.(2.4) makes it possible to test some of the ana- 
lytical limitations of the spiral vortex model, since the two-dimensional enstrophy 

spectrum may be computed from numerical solutions of the two-dimensional Navier- 
Stokes equations rather than depending entirely on asymptotic methods. Further, 
one can use more general two-dimensional solutions. Buntine and Pullin (1988) 
and Pullin and Buntine (1989) used numerical analysis in a similar context. Their 
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approach was somewhat different than that described here. They expanded the 
two-dimensional vorticity in a Fourier series, as in Eq.(2.8), and solved numerically 
for the coefficients u> n (r,t). These functions were then used to evaluate the Bessel 
function integrals required for L:(60), and finally the energy spectrum was obtained 
through the appropriate time integration. The advantage of the present approach is 
flexibility. Having identified the required physical quantity as the two-dimensional 
enstrophy spectrum, it may be computed by any method and processed through 
Eq.(2.4) without further reference to the detailed analytical expansions of Lund- 
gren (1982). 

In the computations reported below, a pseudo- spectral method was employed in 
a vorticity/streamfunction formulation in a square region with periodic boundary 
conditions. Aliasing errors were removed by the 2/3 rule and most of the compu- 
tations were done with 240 2 resolution. The computations were performed at the 
Center for Turbulence Research at NASA Ames Research Center on a Cray YMP. 

The characteristic length in the flows was taken such that the sides of the box 
have dimensionless length equal 2tt units; this makes the wavenumbers be integers. 
The actual unit of length was related to some measure of the initial radius of the 
vorticity distribution, thus the vortex is quite a bit smaller than the box. The 
characteristic velocity was selected such that the dimensionless circulation of the 
vortex was unity. This makes the Reynolds number of the flows be Tjv where T is 
the dimensional circulation of the vortex. The Reynolds number would be about 
an order of magnitude smaller if an average swirling velocity at unit radius were 
used for the characteristic velocity. The dimensionless turn around time is about 
40 units. 

The enstrophy spectrum, given by Eq.(2.5), was approximated by summing all 
values of the squared magnitudes of the discrete Fourier coefficients with y wave- 
numbers in a cylindrical shell of unit wavelength in wave space. This value is then 
assigned to a wavenumber which is the midradius of the shell. 

Two series of computations were done with the spiral vortex solution. In the 
first, the two-dimensional enstrophy spectrum was computed from the analytical 
spiral solution by a fast Fourier transform algorithm. The second computes the 
enstrophy spectrum from a numerical solution of the spectral equations with the 
spiral solution as an initial condition. Comparison thus tests the integrity of the 
spiral solution. 

The spiral solution was taken in the form of a two-sided rollup with two halo 
spirals of amplitude /, with 180° separation, plus an additional central core. The 
specific functions used here are expressed as 


u; 2 (r, 0 , t) = u>o(r) + 2/(r) ^ cos [2 n($ - f2(r)f)] exp (-4n 2 fi'(r) 2 j^ 3 /3) , 


n — i 


h^r/aj 2 exp [-(r/a) 2 ] 


fir) = r 


2(1 + h)na 2 


(3.1) 

(3.2) 


(3-3) 
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(a) (b) 

Figure 1.1. t = 0 , t spira i=50. See caption below. 
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Figure 1.2. <=50 , < sp iral=100. See caption below. 
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Figure 1.3. <=100, < S piral=150. See caption below. 


Q(r) = T 


[l + h - (1 + h + /t(r/«) 2 ) exp (~(r/«) 2 )] 
2(1 + h)na 2 


u ^ [w 0 (r) - 2fi(r)] 

(»') = > 


(3.4) 


(3.5) 


where T is the circulation of the vortex, a is a radial dimension and h is a constant. 
If 0 < h < 1, Q(r) will be monotone decreasing. In dimensionless form T = 1, « = 1, 
v is the reciprocal of the Reynolds number and in these computations h 1/— • 
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FIGURE 1.4, £= 150, £spiral=200. See caption below. 
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FIGURE 1,5. £=200, £» P irai=250. See caption below. 


600 


FIGURE 1.1-1. 5. Analytical spiral solution at selected times. Spiral is the time 
in Eq.(3.1), £ is for comparison with figures (2.1)-(2.5) which have the same initial 
condition. T/u = 25,000. (a) vorticity contours. ( b ) vorticity along a horizontal 
cut through the middle of (a), (c) Enstrophy spectrum versus wavenumber. 



FIGURE 1.6. k 5 / 3 times the three-dimensional energy spectrum calculated by 

using Eq.(2.4) with the analytical spiral solution. 
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FIGURE 2.1. 0. See caption below. 
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FIGURE 2.2. tf=50. See caption below. 
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FIGURE 2.3. t— 100. See caption below. 


The first series, displayed in figures 1.1-1. 6, is computed from the analytical 
solution with T/u — 25, 000. The required discrete Fourier series were computed by 
means of Temperton’s fast Fourier transform algorithms. The figures labeled (a) are 
vorticity contours at 5 times starting with the partially wound state corresponding 
to a dimensionless time Spiral = 50; the figures labeled ( b ) show the vorticity along 
a horizontal cut through the middle of the contour figures. The figures labeled (c) 
are of the enstrophy spectrum. There is a huge peak at low wavenumber which 
is off scale in these figures. Enstrophy similarity, which was found asymptotically, 
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(•) ( 6 ) 

Figure 2.4. 2=150. See caption below. 
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(a) (*) 

FIGURE 2.5. 2=200. See caption below. 



FIGURE 2. 1-2. 5. Computed spiral solution at selected times. T/u = 25,000. 

(a) vorticity contours. ( b ) vorticity along a horizontal cut through the middle of 
(a), (c) Enstrophy spectrum versus wavenumber. 



FIGURE 2.6. fc 5 / 3 times the three-dimensional energy spectrum calculated by 

using Eq.(2.4) with solution computed using the analytical spiral solution as initial 
condition. 
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FIGURE 3.1. <=0. See caption below. 




FIGURE 3.2. t= 50. See caption below. 
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Figure 3.3. f=100. See caption below. 


is not perfect but is qualitatively recognizable. As time increases, one can see the 
number of turns in the spiral increase due to differential rotation and the peak in 
the enstrophy spectrum move outward to larger wavenumber as the spatial scale in 
the spiral decreases. 

Figure 1.6 shows the three-dimensional energy spectrum computed from Eq.(2.4). 
The integration was carried out from T = 0 to T = 250 with a coarse integration 
time interval of AT = 5, which was adequate here. (In computations where the 
Navier-Stokes equations were solved numerically, AT was taken the same as the 
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FIGURE 3.4. £=150. See caption below. 
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FIGURE 3.5. £=200. See caption below. 
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FIGURE 3.6. £=250. See caption below. 


FIGURE 3.1 — 3.6. Computed spiral solution at selected times. T/v = 100,000. 
(a) vorticity contours. ( b ) vorticity along a horizontal cut through the middle of 
(a), (c) Enstrophy spectrum versus wavenumber. 


updating time step, namely AT = .05.) The strain-rate was taken to be unity. The 
low wavenumber end of the spectrum was suppressed by cutting off the integration 
when S(T)~ l / 2 k is less than 5, thus avoiding the large low wavenumber enstrophy 
peak. The result in figure 1.6 has a wavenumber range from about 60 to 120 where 
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k 


FIGURE 3.7. k 5 / 3 times the three dimensional energy spectrum calculated by 

using Eq.(2.4) with the solution depicted in figures 3.1- 3.6. 


the spectrum is approximately k 5 / 3 . Note that since 2 1 / 3 ~ 1,26, one should be 
able to tell the difference between A:"" 5 / 3 and k ~ 2 on this figure. 

The second series of computations, seen in figures 2. 1-2.6, are presented in the 
same format as the first series. The Navier- Stokes equations were solved numerically 
with initial conditions the same as the initial frame of the first series. Here the object 
is to show that the analytical spiral solution is a good approximation to the Navier- 
Stokes equations even though time (f spiral ) is not large enough at the beginning for 
the spiral to have many turns. Comparing the two series of computations, one can 
see that the results are quite close but not identical. In particular the enstrophy 
spectrum of the second series is smoother and more compact. The energy spectrum 
is very similar but perhaps slightly tipped toward fc -2 . 

A third series of computations has been carried out with quite different initial 
conditions with the objective of producing a spiral solution without actually start- 
ing with one. The vorticity distribution at the initial time is shown in figure 3.1. 
One large vortex, with radius 1.5, almost uniform vorticity and circulation .6 is 
surrounded by 10 smaller vortices with centers 2 units from the center of the large 
vortex. Each of these smaller vortices have radii .3 units and circulation .04 so 
that the total circulation of the configuration is unity. The Reynolds number was 
Tjv = 100,000. The smaller vortices get sheared into bands of vorticity which 
continue to tighten in the differential rotation of the combined vorticity. The re- 
sult is a multilayered spiral vortex with a decent enstrophy cascade. Similarity is 
approximately satisfied; comparing the enstrophy spectrum at t = 50 with that at 
t = 100, the amplitude is almost half and the peaks have moved approximately to 
twice the wavenumbers as required by similarity. However, comparison of t = 100 
with t = 200 is not quite as satisfactory. 

The energy spectrum in figure 3.7 has a short Ar~ 5 / 3 range from about k = 36 to 
k — 72 and then falls off faster. The viscous factor in Eq.(2.15), which is about .9 
when k = 120, is not small enough to account for all of the observed decrease. 
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4. Conclusions 

The original 1982 spiral vortex model was reformulated in a form where the effect 
of the two-dimensional flow in the vortex cross section was more clearly separated 
from the axial stretching. While this was done in order to be able to use compu- 
tational methods more effectively, the new formulation has allowed greater insight 
into the workings of the model. It was shown in section 1 that, in the inviscid limit, 
the analytical spiral vortex solution gives an asymptotic time-dependent enstrophy 
spectrum in a self-similar form which conserves enstrophy. It is this self-similar 
form which leads to the A: -5 / 3 energy spectrum. 

Computations with flows which develop spiral structure showed, qualitatively cor- 
rect, but imperfect, self-similar enstrophy spectra. The resulting three-dimensional 
energy spectra nevertheless showed short ranges with the fc -5 / 3 power law. These 
results seem quite rugged and verify results which were previously obtained asymp- 
totically. 

The integral which processes the computed two-dimensional enstrophy spectrum 
to produce the three-dimensional energy spectrum requires integrations over a very 
long times, of the order of 6 turn-around-times, during which the vortex is stretched 
by a factor of about 250. The result from appendix A, which shows that the strain 
may be applied intermittently, make this seem more reasonable. 
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Appendix A. Insensitivity to time dependent strain-rate 

The strain-rate which a real vortex feels is generated by the presence of other 
nearby vortices and is unlikely to be constant for very long; therefore, it is important 
to see if the results are sensitive to time dependence of the strain-rate function. 

Equation (2.4) calls for the stretching function S to be expressed as a function of 
the strained time T. When the strain-rate a(t) is constant this leads to the simple 
linear relationship given by Eq.(2.7). When the strain-rate is not constant, 5 and 
T are defined by 

S(t) = exp ^ j (.41) 

T(t)= f S(t')dt ' (A2) 

Jo 

and it is not simple to relate them analytically, though it is clear that since S is 
positive, T is a strictly increasing function of time which can be inverted in principle. 
It is easy to get the relationship numerically in special cases. As an example, the 
strain-rate function 


a(t) = 1 + .5sin(27r£) 


(.43) 
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FIGURE 4. Time dependent strain-rate, (a), (b) generated with sinusoidal strain- 
rate, a(t) = 1 + .5sin(27r*); (c),(d) generated with on-off strain-rate of same period, 
a(t) = l + (_l) int (20 . 


has average value unity and periodic variations of ±50%. During each successive 
period, the stretch increases by a factor e. These functions have been computed (by 
solving the pair of differential equations dS/dt = aS, dT/dt = S) and presented as 
S vs. T in figure 4a. The relationship is roughly linear. This function was used in 
Eq.(2.4) to compute the energy spectrum using the analytical spiral solution with 
the same set-up a s for the “series one” computations. The resulting energy spectrum 
in figure 4b is almost identical with that in figure 1.6 which had constant strain-rate. 
Therefore, it appears that the spectral results are insensitive to moderate variations 
in strain-rate. 

A more extreme intermittent case was tried in which there are alternating periods 
of positive strain-rate and zero strain-rate. The function used was 

a(t) = l + (— l) int(2<) (.44) 

where the function “int” truncates the decimal part of a number. This strain-rate 
function has the same period as Eq.(A3) and the same average value. Alternating 
half-periods have a = 2 or a = 0. The S vs. T result in figure 4c is still very roughly 
linear, and the energy spectrum still has the same range of fc“ 5 / 3 . 
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As one increases the frequency of the strain-rate function, the relationship be- 
tween 5 and T will become more nearly linear and will have little effect on the 
spectral result. Lower frequencies could have an effect. In the extreme case of very 
long on and off periods, one would get either A:” 5 / 3 or depending on which 
period comes first. 
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